Kinetic energy of protons in ice Ih and water: a path integral study 
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The kinetic energy of H and O nuclei has been studied by path integral molecular dynamics 
simulations of ice Ih and water at ambient pressure. The simulations were performed by using 
the q-TIP4P/F model, a point charge empirical potential that includes molecular flexibility and 
anharmonicity in the OH stretch of the water molecule. Ice Ih was studied in a temperature range 
between 210-290 K, and water between 230-320 K. Simulations of an isolated water molecule were 
performed in the range 210-320 K to estimate the contribution of the intramolecular vibrational 
modes to the kinetic energy. Our results for the proton kinetic energy, Kh, in water and ice Ih show 
both agreement and discrepancies with different published data based on deep inelastic neutron 
scattering experiments. Agreement is found for water at the experimental melting point and in the 
range 290-300 K. Discrepancies arise because data derived from the scattering experiments predict 
in water two maxima of Kh around 270 K and 277 K, and that Kh is lower in ice than in water at 
269 K. As a check of the validity of the employed water potential, we show that our simulations are 
consistent with other experimental thermodynamic properties related to Kh, as the temperature 
dependence of the liquid density, the heat capacity of water and ice at constant pressure, and the 
isotopic shift in the melting temperature of ice upon isotopic substitution of either H or O atoms. 
Moreover, the temperature dependence of Kh predicted by the q-TIP4P/F model for ice Ih is found 
to be in good agreement to results of path integral simulations using ab initio density functional 
theory. 

PACS numbers: 65.20.-w, 65.40.Ba, 61.20.Gy, 82.20.wt 
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I. INTRODUCTION 

The kinetic energy of atomic nuclei in molecules and 
condensed phases in thermodynamic equilibrium is re- 
lated to the interatomic potentials present in the system. 
Thus, although in the classical limit the average kinetic 
energy of a nucleus does not depend on details of its en- 
vironment and interatomic interactions, but only on the 
temperature (equipartition principle), in a quantum ap- 
proach the kinetic energy of a nucleus gives information 
on the effective potential in which the particle moves. 

Deep inelastic neutron scattering (DINS) has been re- 
cently applied to study the thermal average of the kinetic 
energy of protons, Kh, in water and ice. Kh is derived 
from the measured neutron scattering function by data 
analysis that includes several numerical corrections, as 
the subtraction of the cell signal, 1 the application of the- 
oretical models as the impulse approximation (that con- 
siders that the neutrons are scattered by single atoms 
with conservation of momentum and kinetic energy of 
the neutron plus the target atom), and numerical proce- 
dures to perform an integral transform of the scattering 
function to obtain the momentum distribution of the tar- 
get atomi2 Therefore analysis of experimental DINS data 
is not free of numerical limitations that have led in the 
past to reinterpretation of results. For example, the ef- 
fective potential of H along the stretching direction of 
the OH bond in supercritical water at 673 K was con- 
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sidered as a double well potential on the basis of DINS 
measurements, 3 however later experiments with longer 
counting times and subtraction of the cell signal did not 
show any evidence for such double well potential. 4 

For water, DINS investigations have been performed at 
atmospheric pressure and various temperatures: 300 K, 4 
296 271 K£ and 269 K£ The last two temperatures 
correspond to supercooled water and a substantial in- 
crease in K H of about 58% (at 271 K) and 40% (at 269 K) 
was reported with respect to ambient water. This large 
increase was heuristically interpreted to be a consequence 
of a double well potential felt by delocalized H atoms. 5 - 5 
However, a controversy has arisen about the data analy- 
sis at these temperatures, and it has been claimed that 
the observed data can be also explained without invoking 
large changes in Kh and without assuming the dereal- 
ization of H in a double well potential^ Very recently 
a new set of DINS experiments have been performed for 
water at temperatures in the range 272-285 K, and seem- 
ingly what appears as a unique set of measurements have 
been published by two different groups. 8,9 From an analy- 
sis of their experimental data these authors found a max- 
imum in Kh at about 277 K that was considered as an 
evidence of a new water structural anomaly related to 
the well-known existence of a maximum in the density of 
water at the same temperature. 8 DINS derived values of 
Kh have been also reported at 5 K and at 269 K for the 
stable phase of ice at atmospheric pressure (ice Ih)i2^£ 

The calculation of Kh in ice Ih has been done so far 
with the help of empirical models consisting of a set of 
decoupled quantum harmonic oscillators whose frequency 
is derived from optical data and measured vibrational 
density of states. The contribution of each harmonic os- 
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cillator to K h was determined by the average number of 
excited phonons at the given temperature and by an es- 
timation of the H participation in each normal mode.-U> 
This numerical approach predicts that Kh in ice Ih is 
nearly constant between 5 and 300 K. Similar empirical 
calculations have been performed for water^ A limita- 
tion of such approaches is that they are based on ex- 
perimental information and neither the pressure nor the 
density of the condensed phase can be explicitly consid- 
ered in the model. 

Momentum distributions of protons in ice and water 
have been derived by path integral (PI) simulations spe- 
cially designed to sample non-diagonal density matrix el- 
ements. The first attempts used empirical water models 
to study ice at 269 K t 12 i 13 and supercritical water at 673 
KM- Water molecules were treated later by ab initio den- 
sity functional theory (DFT) in studies of ice at 269 K 
and water at 300 K^i obtaining for the radial proton 
momentum distribution closer agreement to experiment. 
The thermal average of the proton kinetic energy Kh, 
can be calculated from the width of its momentum dis- 
tribution. However, Kh is more easily accessible in stan- 
dard PI simulations, that sample diagonal elements of the 
density matrix, by the calculation of the virial estimator 
of the kinetic energy 15 i 16 PI simulations of the kinetic 
energy of atomic nuclei have shown good agreement to 
data derived from neutron scattering experiments in solid 
neoni 17 ' 18 

In the present paper the virial estimator of the kinetic 
energy, K, has been calculated for water and ice Ih at 
atmospheric pressure by PI simulations. The empirical 
q-TIP4FP/F model has been chosen for the water sim- 
ulations because it is an anharmonic flexible potential 
that reproduces the experimental density-temperature 
curve of water at atmospheric pressure with reasonable 
accuracy^ This potential model has been recently em- 
ployed to study the isotopic shift in the melting temper- 
ature of ice ihi 2 ^ The partition of the kinetic energy into 
H- and O-atom contributions (Kh, Ko) has allowed the 
comparison to data derived from DINS experiments. To 
assess the reliability of our simulations for Kh we dis- 
cuss also several thermodynamic properties that depend 
on the kinetic energy, as the constant pressure heat ca- 
pacity, C p , of ice and water and isotopic shifts in the 
melting temperate of ice. An additional check of the 
employed potential model has been conducted by com- 
paring the temperature dependence of Kh in ice Ih to 
results derived by PI simulations based on ab initio DFT 
calculations using the SIESTA methodj 21 ' 22 

The structure of this paper is as follows. In Sec. |TT] 
a short summary of the employed computational condi- 
tions is presented. Then the results obtained with the 
q-TIP4P/F model for the thermal average of the kinetic 
energy of ice, water, and an isolated molecule are given 
and compared to available experimental data and also to 
ab initio DFT simulations. In Sec. Mil the simulation 
results of C p at 271 K are presented for ice and water 
and compared to experiment. In Sec. IIVI it is shown 



that kinetic energy differences between solid and liquid 
phases at isothermal conditions determine the sign of the 
isotopic shift in the melting temperature. Finally, we 
summarize our conclusions in Sec. [V] 



II. KINETIC ENERGY 
A. Computational conditions 

In the PI formulation of statistical mechanics the parti- 
tion function is calculated through a discretization of the 
integral representing the density matrix. This discretiza- 
tion defines cyclic paths composed by a finite number L 
of steps, that in the numerical simulation translates into 
the appearance of L replicas (or beads) of each quantum 
particle. Then, the implementation of PI simulations 
relies on an isomorphism between the quantum system 
and a classical one, derived from the former by replac- 
ing each quantum particle (here, atomic nucleus of H 
and O atoms) by a ring polymer of L classical particles, 
connected by harmonic springs with a temperature- and 
mass-dependent force constant. Details on this compu- 
tational method are given elsewhere i 2 ^— The configu- 
ration space of the classical isomorph can be sampled 
by a molecular dynamics (MD) algorithm, that has the 
advantage against a Monte Carlo method of being more 
easily parallelizable, an important fact for efficient use of 
modern computer architectures. Effective reversible in- 
tegrator algorithms to perform PIMD simulations have 
been described in detail in Refs. [27i - [30l Ref. [U in- 
troduces useful algorithms to treat full cell fluctuations 
and multiple time step integration. All calculations were 
done using originally developed software and paralleliza- 
tion was implemented by the MPI library^ 

PIMD simulations in the isothermal-isobaric NPT en- 
semble (N being the number of particles, P the pres- 
sure, and T the temperature) were conducted for ice Ih 
and water by using the point charge, flexible q-TIP4P /F 
model. This model was parameterized to provide the 
correct liquid structure, diffusion coefficients, and in- 
frared absorption frequencies (including the translational 
and librational regions) in quantum simulations. 19 Wa- 
ter simulations were done on cubic cells containing 300 
molecules, while ice simulations included 288 molecules 
in a proton disordered orthorhombic simulation cell with 
parameters (4a, 3-\/3a, 3c), with (a, c) being the standard 
hexagonal lattice parameters of ice Ih. The total kinetic 
energy is defined as K = 2Kh + Ko- The kinetic energy 
of the H and O nuclei (Kh, Kq) were derived by the 
virial estimator ^ 5 ' 16 ' 29 We stress that Kh and Ko are 
obtained as expected averages using the exact quantum 
partition function of the q-TIP4P/F water model, i.e., 
our calculation of quantum kinetic energies is free of any 
other physical assumption apart from the statistical un- 
certainty inherent to the actual PIMD simulations and 
the use of an empirical potential model. Expected aver- 
ages were derived in runs of 5 x 10 5 MD steps (MDS) for 
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water and 2.5 x 10 5 MDS for ice Ih, using in both cases 
a time step of 0.3 fs. The presence of disorder and spa- 
tial density fluctuations in the liquid causes that larger 
simulation runs are required for water. The system equi- 
libration was conducted in runs of 5 x 10 4 MDS. To have 
a nearly constant precision in the PI results at different 
temperatures, the number of beads L was set as the in- 
teger number closest to fulfill the relation LT — 6000 K, 
i.e., at 210 K the number of beads was L = 28. Ad- 
ditional computational conditions are identical to those 
employed in Ref. [2^ and they are not repeated here. The 
simulations were done at atmospheric pressure in a tem- 
perature range of 210-290 K for ice Ih and 230-320 K for 
water. 

Additionally the kinetic energy, Ki ntra , associated to 
the three vibrational modes of an isolated water molecule 
was obtained by PIMD simulations at temperatures be- 
tween 210-320 K. In the case of an isolated molecule the 
virial estimator allows us the calculation of the oxygen, 
K ,intra, and hydrogen, K Htintra , contributions to the 
vibrational energy, Kintra- The six translational and ro- 
tational degrees of freedom of an isolated water molecule 
behave classically at the studied temperatures and their 
kinetic energy (fcsT/2 per degree of freedom) has to be 
summed to Ki ntra to obtain the total kinetic energy, K, 
of the molecule. If mo and M are the masses of the 
oxygen atom and the water molecule, then the fraction 
of translational kinetic energy corresponding to the oxy- 
gen atom is given by mo/M. Analogously for a classical 
rotation around a given axis, the fraction of the kinetic 
energy corresponding to the O atom is given by Io / I, 
where Io and / are the moment of inertia of the O atom 
and the water molecule with respect to the rotation axis. 
This ratio has to be calculated for three principal axes 
of the molecule to obtain the partition of the rotational 
kinetic energy into O and H contributions. 33 The energy 
unit used throughout this work (k J / mol) refers to a mole 
of either the molecule (H2O) or atom (O or H) under 
consideration. 

Finally, we have performed PIMD simulations in ice Ih 
by coupling our PIMD program to the SIESTA code ) 21 i 22 
in such a way that the potential energy, atomic forces, 
and stress tensor used in the NPT simulations are de- 
rived from ab initio DFT calculations. The present im- 
plementation allows us to study only small simulation 
cells, that are not large enough for an appropriate de- 
scription of the disorder in the liquid state. We ex- 
pect to overcome this limitation in the future by al- 
lowing the simultaneous parallelization of both PI and 
electronic structure codes. DFT results were derived for 
a 2 x 1 x 1 supercell of the standard hexagonal cell of 
ice Ih with the aim of comparison to the data obtained 
with the q-TIP4P/F model. The DFT calculations were 
done within the generalized-gradient approximation with 
the PBE functional. 34 Core electrons were replaced by 
norm conserving pseudopotentials 35 in their fully non- 
local representation. 36 In this study a double-^ polarized 
basis set was used and a grid of 27 k-points was employed 
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Figure 1: Thermal average of the kinetic energy of ice Ih and 
water as a function of temperature, as calculated by PIMD 
simulations with the q-TIP4P/F potential model at P = 1 
atm. The kinetic energy, K, of an isolated molecule and its vi- 
brational contribution, Ki ntr a, are also given. The calculated 
error bars for ice and water are smaller than 0.01 kj/mol and 
for the molecule amount to 0.05 kj/mol. 



for the sampling of the Brillouin zone of the solid. In- 
tegrals beyond two-body interactions are performed in a 
discretized real-space grid, its fitness determined by an 
energy cutoff of 100 Ry. Average properties were derived 
by PIMD simulations runs of 4 x 10 4 MDS using a time 
step of 0.5 fs. Classical MD simulations of water using 
the SIESTA code have been published in Refs. UTtM 



B. Total kinetic energy 

The temperature dependence of K calculated for ice 
Ih, water, and an isolated molecule with the q-TIP4P/F 
model potential is represented in Fig. [T] We observe 
that at a given temperature K is always larger in ice Ih 
than in water. For a molecule K is the sum of global rota- 
tional and translational terms (which are classical contri- 
butions at the studied temperatures), and a vibrational 
term, Ki ntra , corresponding to the two OH stretches and 
the HOH bending mode. In Fig. Q] the temperature de- 
pendence of the Ki n tra term has been also shown for the 
isolated molecule. While K in ice and water increases 
with temperature, the vibrational term, Ki ntra , of the 
molecule remains nearly constant. This fact indicates 
that the three intramolecular vibrational modes remain 
mainly in their ground state in the studied temperature 
range, while vibrational modes related to the H-bond net- 
work display some amount of thermal population of ex- 
cited states. 

The partition of K in ice Ih and water into intra- and 
intermolecular contributions (Kintra and K inter) is of in- 
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terest to quantify the effect of the presence of H-bonds 
in this observable. The intramolecular term, K intra , 
is attributed to the internal vibrational modes of the 
molecule, while the intermolecular one, Ki nter , is related 
to the librational (rotational) and translational motions 
of the water molecules in the condensed phases. The 
presence of a H-bond network causes a large broadening 
of the internal vibrational modes of the water molecule 
and a shift in the frequencies of their peaks in the in- 
frared absorption spectra. 19 This fact implies that both 
inter- and intramolecular motions are clearly coupled in 
the condensed phases of water. Therefore the partition of 
the total kinetic energy is presented here only as an ap- 
proximation for didactic purposes in ice and water, and 
it is not intended to suggest that both types of degrees 
of freedom may be treated as separable variables. 

We estimate Ki ntra in ice Ih and water with the help 
of a proportionality relation based on the fact that at the 
studied temperatures Ki ntra is close to the kinetic energy 
of the ground state associated to three intramolecular 
modes, Kcs.intra- Particularizing to the case of ice we 
have, 

Kintrajice) ^ K GStintra {ice) 
K mtra (molecule) K GS ^ n tra{molecule) 

As a rough estimation of the ground state kinetic energy, 
KGs,intra, in either ice, water or an isolated molecule, 
we employ here a quasi-harmonic approximation, so that 
a vibrational mode of wavenumber uj displays a kinetic 
energy of hw/A in its ground state. Thus Kcs,intra will 
be approximated as 



Kos,intra ~ J {^■ u3 OH + ^HOHJ ■ (2) 

The wavenumber oj gh is used here as an estimation of the 
average wavenumber of the two stretching modes (sym- 
metric and asymmetric) of an isolated water molecule 
or the corresponding vibrational bands in the con- 
densed phases, is determined from the equilibrium 
bond distance, dou, of either ice, water or an isolated 
molecule, by calculating the second derivative of the po- 
tential energy with respect to doH, and by considering 
the actual O and H masses. 20 The bending wavenum- 
ber is ujhoh =1600 cm" 1 . We recall that the bending 
potential is described by a simple harmonic term in the 
employed water model. 19 The ground state assumption 
for Ki n tra in Eq. (fTJ) is further justified by noting that the 
excitation energy corresponding to the vibrational mode 
of lowest frequency, Hujhoh-, is about seven times larger 
than the thermal energy, ksT, at the highest studied 
temperature (320 K), and thus the thermal population 
of excited vibrational modes, as given by a Bose-Einstein 
statisticS) 33 i 40 results vanishingly small. 

The partition of K into Ki ntra and Ki nter contribu- 
tions is presented in Table U at a temperature of 270 K. 
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Figure 2: Oxygen contribution, Ko, to the kinetic energy of 
ice Ih and water as a function of temperature. The classical 
limit for ice and water amounts to (3/2)fcsT per atom. The 
calculated error bars are smaller than 3xl0 -3 kj/mol. For a 
given temperature, the difference of Ko in ice Ih and water 
(~ 0.05 kj/mol) is significantly larger than the error bar. 

For water the intermolecular contribution, Ki nteri is es- 
timated to be 30% (10.1 kJ/mol) of K, while for ice Ih 
this contribution is slightly larger (32%, 10.7 kJ/mol), 
in agreement with the expectation that H-bonds in ice 
are stronger than in water 41 We display also the average 
OH bond distance, doHi an d the corresponding quasi- 
harmonic stretching wavenumber, w^. The bond dis- 
tance increases in the order molecule <§C water < ice. We 
note that x-ray Compton scattering studies of water and 
ice Ih have shown that the elongation of the OH bond in 
ice Ih is correlated to its stronger H-bond network^ 2 - in 
agreement to the simulation results for don and Ki nter 
shown in Tab. UJ Our result of K in ice Ih displays a 
reasonable agreement to a previous calculation that is 
based on the use of experimental vibrational frequencies 
assuming a harmonic approximation (see the last column 
of Tab. [E])4I 



C. Kinetic energy of O and H nuclei 

The O nuclei contribution, Ko, to the kinetic energy 
is presented in Fig. [5] Ko in the condensed phases is 
significantly larger than its classical limit, i.e., quantum 
effects related to the oxygen mass are relevant for its 
kinetic energy. At a given temperature Ko is slightly 
larger in ice Ih than in water (see also the Ko values of 
ice and water at 270 K in Table fl}. 

The temperature dependence of Kh is presented in 
Fig. [3] for ice and water. The Kh contribution to the 
kinetic energy is substantially larger than the Ko term 
and it is also much larger than the classical limit, i.e., 
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Table I: Thermal average of the kinetic energy, K, of water and ice Ih at atmospheric pressure and of an isolated water molecule. 
The results were derived by PIMD simulations at 270 K. The energy K has been partitioned into O- and H-atom contributions 
(Ko, Kh)- For water and ice the estimated intramolecular and intermolecular contributions to K, Ko, and Kh are also given. 
For a molecule the classical translational and rotational energy is given as K tr + r ot- This classical term has been partitioned 
into O and H components, don is the thermal average of the intramolecular OH distance, and ijJq H is the quasi-harmonic 
stretching wavenumber calculated at the distance don- The last column displays the theoretical results of Ref. [ll] for ice Ih at 
269 K. Kinetic energies are given in kj/mol (to convert to meV multiply by a factor of 10.37). The standard error in the final 
digit is given in parenthesis. 
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Figure 3: Proton contribution, Kh, to the kinetic energy of 
ice Ih and water as a function of temperature. The classical 
limit for the condensed phases corresponds to (3/2)fesT per 
atom. The calculated error bars are smaller than 5 x 10~ 3 
kj/mol. At constant temperature, the difference of Kh in ice 
Ih and water (~ 0.1 kj/mol) is significantly larger than the 
error bar. 



quantum effects related to their low atomic mass are very 
important for H atoms. At a given temperature Kh is 
always larger in ice Ih than in water. 



The partition of Kh and Ko into intramolecular and 
intermolecular contributions can be derived by assuming 
the following proportionality relation for the intramolec- 
ular kinetic energy in the condensed phases 



K 



H. intra 



Ko 



intra 



K. H,intra 
Ko .intra 



(3) 



molecule 



i.e., the ratio K h, intra / Ko^ntra at a given temperature 
is assumed to be constant (~ 7.7) for ice, water, and 
an isolated molecule. Although we can not provide a 
theoretical estimation of the error of this assumption, 
water molecules remain as clearly recognizable entities 
in the condensed phases at atmospheric pressure, and 
therefore the intermolecular interactions are not expected 
to change drastically the constant ratio assumed by Eq. 
^ . In Table |T] the result of this partition is given at 
270 K. For water, the intermolecular contribution to Ko 
is larger than the intramolecular term (71% and 29%, 
respectively). However in the case of Kh in water, the 
intermolecular contribution (23%, 3.3 kj/mol) is much 
lower than the intramolecular term (77%, 11 kj/mol, cor- 
responding to the two OH stretches and the HOH bend- 
ing mode). For ice Ih, having a stronger H-bond network 
than water^i we find that intermolecular contributions 
to Ko and Kh are slightly larger than in water, while 
the intramolecular contributions are slightly lower. 
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Figure 4: PIMD results obtained with the q-TIP4P/F model 
for the proton kinetic energy in water and ice Ih are compared 
to data derived from DINS experiments at temperatures in 
the range 269-300 Ki2A2riS Also PIMD results for ice Ih de- 
rived by ab initio DFT calculations with the SIESTA code 
are given. Lines are guides to the eyes. 



D. Comparison to DINS data 

In Fig. [J] our Kh results are compared to data derived 
from DINS experiments^— Our calculation for water 
agrees reasonably well with the data derived from the 
neutron measurements at 273 K and the three values in 
the range 290-300 K. However, significant deviations are 
found in the temperature ranges between 269-272 K and 
274-286 K. 

The value of K H derived at 271 K from the DINS ex- 
periment predicts an excess of 58% (8 kJ/mol) with re- 
spect to the room temperature result. This excess was 
associated to a coherent derealization of the proton be- 
tween neighbor oxygen atoms. 8,9 However, this tentative 
explanation has been criticized by noting that measured 
proton momentum distribution implies a contracted pro- 
ton wave function which is inconsistent with the pur- 
ported increased proton derealization^ Moreover, co- 
herent derealization (i.e., tunneling) of a proton is as- 
sociated to a resonance of the particle between two lo- 
calized states. A coherent derealization of the proton 
at T c —271 K would imply that the resonance frequency 
associated to the coherent process must be large enough 
to display the following relation to the thermal energy: 
Huj c 3> ksT c . In other case the quantum coherence would 
be destroyed by thermal excitations. This inequality 
should remain valid if T c increases to 273 K, i.e., a slight 
raise of 2 K should not provide enough thermal energy 
to destroy a coherent process of wavenumber uj c for the 
proton. Therefore, it seems somewhat contradictory that 
the excess of Kh attributed to a coherent process at 271 
K disappears completely at 273 K. 



An alternative explanation for the large value of Kh 
as a localized proton state at 271 K leads also to an un- 
physical picture. In this respect, as long as the water 
molecules maintain their molecular structure, i.e., with- 
out a dramatic increase in their stretching frequencies 
(which has been never observed) , the intramolecular con- 
tribution to Kh must remain below the isolated molecule 
result of about 11 kJ/mol (see Fig. [3]). Therefore an ex- 
cess as large as 8 kJ/mol in the intramolecular energy 
seems to be incompatible with the experimental molecu- 
lar structure of water. 

The second maxima in K h derived from the DINS ex- 
periments in the temperature range 274-286 K has been 
attributed to a new structural anomaly of water related 
to the existence of a maximum in the water density as 
a function of temperature^^ This explanation is how- 
ever not supported by our simulations. The employed 
water potential is able to reproduce realistically the ex- 
perimental density maximum of water as a function of 
temperature (see Fig. [S]) ^ In the temperature range be- 
tween 274-286 K there appears a maximum in the density 
of water. Both experimental and simulation data agree 
in the fact that the water density in this temperature 
range changes by less than 1/1000. Our simulations pre- 
dict that such a small density change does not affect K h 
in any significant way (see Fig. 0}. However, the excess 
in Kh at 277 K as derived from the DINS experiments 
amounts to 2.2 kJ/mol with respect to the room temper- 
ature value. This excess energy is rather large, and of 
similar magnitude to the total H-bond network contribu- 
tion to Kh, that amounts to about 3.3 kJ/mol for the 
employed potential model (see Table U for our result of 
KH,inter hi water at 270 K, this value remains nearly con- 
stant in the studied temperature range). Thus an excess 
of 2.2 kJ/mol represents, for the employed water model, 
an increase of 66% in the kinetic energy associated to the 
intermolecular interaction. We believe that such a large 
energetic modification should imply a dramatic change of 
the H-bond network that however has not been observed 
in any other spectroscopic or thermodynamical data of 
water around 277 K. 

Another discrepancy between our simulations and the 
DINS analysis of Kh is that we find that at isothermal 
conditions Kh is about 0.14 kJ/mol larger in ice than in 
water. However the DINS value reported for ic&iS at 269 
K is much lower (~ 7 kJ/ mol) than that of water— at the 
same temperature. 

It is difficult to find a definite explanation for the agree- 
ment found between simulation data for Kh in water and 
those derived from neutron scattering at some tempera- 
tures (273 K and 290-300 K) and the discrepancies found 
at other temperatures (269-272 K and 274-286 K) and 
in the result for ice Ih. Both theoretical and experimen- 
tal approaches include basic assumptions that might lead 
to unphysical results. From a theoretical side, the main 
limitation comes from the use of an empirical water po- 
tential. In particular, any dynamic process involving de- 
localization of protons between different molecules (i.e., 
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Figure 5: Density of water at 1 atm obtained from PIMD 
simulations and experimental data from Ref. |46|. Simulation 
results correspond to runs of 4 x 10 6 MDS. The continuous 
line is a cubic fit to the simulation data. 



breaking and forming of intramolecular OH bonds) can 
not be described by the employed q-TIP4P/F potential. 
From the experimental side, the bibliography shows an 
experimental controversy concerning the validity of the 
models employed for the data analysis of the DINS ex- 
periments at 269 and 271 Kj&£ Moreover, the numerical 
data analysis of the diffraction experiment implies correc- 
tions due to the experimental setup (detector position, 
subtraction of the cell signal) and a numerical integral 
transform that might be subject to numerical impreci- 
sions. It has been mentioned in the Introduction that 
numerical reinterpretation of DINS experiments have al- 
ready happened in recent timesj^ 

In the face of this disagreement, the rest of the pa- 
per is intended as a consistency check for the ability of 
the employed empirical potential to reproduce the kinetic 
energy of the atoms in ice Ih and water. Firstly, we go 
beyond the empirical potential by presenting results of 
ab initio DFT simulations. 



E. Beyond the empirical potential 

In Fig. @] we show the Kh results for ice Ih derived 
from PIMD simulations using atomic forces calculated 
with the SIESTA codei 21 ' 22 The temperature dependence 
of Kh predicted by the q-TIP4P/F model shows good 
agreement to the ab initio DFT results. The main dif- 
ference between both sets of results is that values pre- 
dicted by the q-TIP4P/F model are ~ 3% larger than 
the DFT ones. The origin of this shift is related to a dif- 
ferent description of the OH stretching, as evidenced by a 
normal mode analysis that provides harmonic stretching 
frequencies about 4% larger in the case of the q-TIP4P/F 



Table II: Experimental results for C p of water (Ref. |44T ) and 
ice Ih (Ref. l4a ) are compared with PIMD results at 271 K and 
1 atm. The calculated C p has been decomposed into potential 
(C p i and C V 2) and kinetic (C v z and C v a) energy contributions. 
The numerical estimation of C P 4 from the data derived from 
DINS experiments^^ is also given. The employed unit is 
Jmol^KT 1 . 





water 


ice 


C p experimental 


76.1 


37.4 


C p PIMD 


71.4±2 


36.1±0.8 




-9.0±0.4 


-7.4±0.2 


Cp2 (Vinter) 


66.3±1 


31.6±0.5 


C p3 (Ko) 


10.3±0.2 


9.9±0.05 


Cp 4 (K H ) 


3.7±0.2 


2.0±0.05 


Cp 4 (K H , DINS) 


-10 4 





model. The ab initio data were derived in a small sim- 
ulation cell (see Sec. Ill A)) . Thus, the finite size error 
has been estimated by comparing Kjj values obtained 
in a small (2a, a, c) and a large (4a, 3\/3a, 3c) supercell 
(i.e., containing 8 and 288 water molecules, respectively) 
by using the q-TIP4P/F model. This finite size effect 
is rather low, the smaller cell shows a Kh value 0.3% 
lower than the larger one. We stress that the q-TIP4P/F 
potential was parameterized to give correct structure, dif- 
fusion coefficient, and infrared absorption frequencies in 
the liquid stated 9 Therefore the reasonable agreement to 
the ab initio data found for Kh in ice Ih should be con- 
sidered as a realistic prediction of the empirical potential 
for a property that was not implicitly built in by the 
potential parameterization. 

In the following Sections we continue our check of the 
q-TIP4P/F potential by comparing some thermodynam- 
ical properties (heat capacity and isotopic shifts) that 
depend on the kinetic energy of the nuclei with experi- 
mental data. 



III. HEAT CAPACITY AT 271 K 

The enthalpy, H, calculated with the employed water 
model can be expressed as 



H ViTjtra ~\~ 'inter 

+ K + K H + PV , 



(4) 



where the first four summands give the internal energy 
as a sum of potential (V) and kinetic energy (K) terms 
and the potential energy is partitioned into a sum of in- 
tramolecular (Vi n tra) and intermolecular (Vi n ter) contri- 
butions. We have fitted our PIMD results for the en- 
thalpy of ice and water to a quadratic function of temper- 
ature in the range 250-290 K. Then the constant pressure 
heat capacity 



r - I — ^ 

p ~\dTJ 



(5) 
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has been calculated at 271 K by performing the tem- 
perature derivative of the fitted function. The results 
are presented in Table [TT| and compared to experimental 
data. The employed water model provides realistic val- 
ues of C p for both water and ice at 271 K. The deviation 
with respect to experiment amounts to about 6% in wa- 
ter and to 4% in ice Ih, with simulated values lower than 
experimental ones. Note that the experimental heat ca- 
pacity of ice Ih is nearly one half of that of water, a fact 
that is reasonably reproduced by the simulations. This 
overall agreement encourages us to analyze separately the 
potential (C P \,C P 2) and kinetic (G P 3,G P 4) energy contri- 
butions to C p , as derived from the temperature derivative 
of the first four summands in Eq. (j4} . At atmospheric 
pressure the temperature derivative of the PV term gives 
a vanishingly small contribution to C p . 

The four contributions to the calculated C p have been 
summarized in Table [XT] We note that the contribution 
C p i of the intramolecular potential energy (Vi n t ra ) is neg- 
ative for both water and ice Ih. This result is an anhar- 
monic effect related to the reduction of the don distance 
with increasing temperature, a fact that has been exper- 
imentally reported in water t 42 i 43 and that is reproduced 
by the employed model potential. The largest contribu- 
tion {C P 2) to the heat capacity comes from the potential 
energy of the H-bond network, Vinter- C P 2 amounts to 
93% of the calculated C p in water and to 88% in ice, a re- 
sult in line to the expectation that the high heat capacity 
of water is a consequence of its H-bond network. Kinetic 
energy contributions (C P 3,Cpi) are comparatively much 
lower. C P 3, derived from the temperature derivative of 
the kinetic energy Ko, is several times larger than C P 4, 
that is derived from Kh ■ In particular, we find that G P 4 
amounts to less than 4 J mol — 1 K _1 (5% of the calculated 
C p ) in water. 

The last line of Table HQ summarizes the value of G P 4, 
as derived from the temperature derivative of the Kh 
values obtained from the DINS experiments, that were 
plotted in Fig. [U This contribution is negative and its 
absolute value (~ 10 4 J mol _1 K _1 ) is several orders of 
magnitude larger than the experimental heat capacity. 
Both sign and magnitude of this estimation of G P 4 seem 
to be unphysical and incompatible with the experimental 
C p data at 271 K. This fact let us to question if the 
maximum in Kh derived from the DINS analysis around 
270 K (see Fig. Ql corresponds to any real physical effect. 



IV. ISOTOPIC SHIFT IN THE MELTING 
TEMPERATURE 




temperature 

Figure 6: Schematic representation of the free energy differ- 
ence, AG, between solid and liquid phases of water as a func- 
tion of temperature at constant pressure. The continuous line 
represents AGh for phases containing the isotope 1 H and the 
melting temperature of the solid is Th- Upon isotopic substi- 
tution of H by deuterium, the free energy curve is shifted to 
the broken one, AGd- At given temperature the shift of AGd 
with respect to AGh (shown by an arrow at the temperature 
Th) can be calculated by the definite integral in Eq. l[7|. A 
normal isotopic shift implies that the melting temperature of 
the heavier isotope is higher (To > Th)- 



ature are a consequence of the dependence of the Gibbs 
free energy, G, on the isotope mass. At constant (T,P), 
this dependence results to be a function of the kinetic en- 
ergy of the isotope. Thus if Gh is the Gibbs free energy 
of a water phase containing 1 H atoms, the free energy of 
the corresponding deuterated phase, Gd, is given by the 
following relation^ 



G 



D 



G 



H 



dm 



K H (m) , 



(6) 



where ttlh and mp are the H and deuterium (D) masses, 
and Kh (m) is the thermal average of the kinetic energy of 
a H isotope with mass m. Note that the mass m appears 
as an auxiliary variable in Eq. ([6]) and corresponds to 
a real isotope only at the integration limits (mff,mi)). 
The Gibbs free energy difference, AG = G s — G/, between 
the solid (s) and liquid (I) phases at constant (T, P) can 
be expressed from Eq. as 



r m ° dm 

AGd - AG ff - / —AK H (m) , (7) 
Jm H m 

where 



For ice Ih at atmospheric pressure experimental iso- 
topic shifts indicate that heavier isotopes display higher 
melting temperatures. Thus deuterated ice Ih melts at a 
temperature 3.8 K higher than normal ice, while the iso- 
topic shift for tritiated ice amounts to 4.5 K. Also H2 18 
melts at a temperature 0.3 K higher than ice Ih made 
of normal H2 16 0. Isotopic shifts in the melting temper- 



AK H (m) = K H , s {m) - K H ,i{m) ■ ( 8 ) 

Kn.si'm) is the isotope kinetic energy in the solid and 
KH.i{fn) the corresponding quantity in the liquid phase. 
A schematic representation of the free energy difference 
AGh at constant pressure is shown in Fig. [B]as a function 
of T. At the melting temperature, Th, the free energy 
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of both phases must be identical and AGh vanishes. At 
temperatures below Tjj we find that AGh < 0, i.e., the 
solid is the thermodynamically stable phase as it displays 
lower free energy than the liquid. However, for T > Th 
one finds AGh > and then the stable phase is the liq- 
uid one. The isotope effect in the melting temperature 
can be obtained by calculating AGd with the help of 
Eq. ||7J| . A normal isotope effect implies that the definite 
integral in Eq. takes a positive value. Thus at con- 
stant temperature AGd becomes smaller than AGh ob- 
taining the schematic result represented in Fig. [51 Note 
that the melting point for the heavier isotope, Tjj (where 
AGd = 0), is shifted towards higher temperatures. 

We conclude that it is the sign associated to the value 
of the definite integral in Eq. ([7]) what determines if the 
isotope effect in the melting temperature is normal (i.e., 
higher melting temperature for a heavier isotope mass 
if the sign is positive) or it is an inverse isotope effect 
(lower melting point for a heavier isotope mass if the sign 
is negative). A sufficient (but not necessary) condition 
to obtain a normal isotope effect is that 

Aifff(m) > 0; for m# < m < mjj (9) 

i.e., the isotope kinetic energy is larger in the solid than 
in the liquid phase for the isotope masses m in the studied 
interval. 

Our simulations show that the condition in Eq. @ is 
satisfied for the isotopic substitution of H by either deu- 
terium or tritium. The analogous condition, AKo(jn) > 
0, is also satisfied for the isotopic substitution of 16 O by 
18 0, with m( 16 0) < m < m( 18 0). Thus the employed 
water model predicts normal isotopic shifts in the melt- 
ing temperature upon isotopic substitution of either H 
or 16 by a heavier isotope, in agreement to experiment. 
Detailed results for deuterated and tritiated phases have 
been presented elsewhere. 2 - The values collected in Ta- 
ble U at 270 K show that for hydrogen we get the value 
AKh(itih) — 0.14 kJ/mol, i.e., Kh is larger in the solid 
than in the liquid, while for 16 we find than Ko is 0.05 
kJ/mol larger in ice than in water. 

The Kh results derived from DINS data for water and 
ice Ih at 269 K give AKh{itih) — —6.7 kJ/mol, i.e., Kh 
is much lower in ice than in water. Note that here the 
inequality of Eq. © is not satisfied for mjj. This fact 
implies that a small increase in the mass of hydrogen 
would produce an inverse isotope effect in the melting 
temperature. This behavior for ran does not imply an 
inverse isotope effect for deuterium, because the latter ef- 
fect is determined by the definite integral in Eq. t na t 
involves all masses m with win < m < mp- However, 
the value of AKh(itih) — — 6.7 kJ/mol at 269 K derived 
from the DINS analysis contrasts in sign and magnitude 
with our result of Aifff (m#) = 0.14 kJ/mol. We stress 
that the employed potential model provides a realistic 
prediction of the isotopic shift in the melting tempera- 
ture of ice, a fact that let us expect that our calculation 
of AKnimn) must be also reasonably realistic^ This 



consideration rises further doubts about the physical sig- 
nificance of the maximum in K h derived from the DINS 
analysis for liquid water around 270 K. 



V. CONCLUSIONS 

Quantum PIMD simulations at atmospheric pressure 
show that the q-TIP4P/F model provides realistic results 
for several thermodynamic properties of water and ice Ih. 
In particular, satisfactory agreement to experiment has 
been found in the dependence of the water density with 
temperature^ the isotopic shifts in the melting temper- 
ature of deuterated and tritiated ice^H and the constant- 
pressure heat capacity of water and ice Ih at 271 K. In the 
present paper we have focused on the calculation of the 
proton kinetic energy, Kh, and its comparison to avail- 
able data derived from the analysis of neutron scattering 
experiments pT— &r— and also to results based on ab initio 
DFT simulations. 

This comparison has shown a remarkable agreement 
between the computational model and DINS data for wa- 
ter at some temperatures (273 K and in the range 290- 
300 K), but a strong disagreement at other similar tem- 
peratures (in the range 269-272 K and 274-286 K). Our 
simulations predict that Kh remains nearly constant in 
the temperature range where DINS data are available at 
atmospheric pressure. On the contrary, Kh values de- 
rived from the scattering experiments display a remark- 
able temperature dependence with two maxima at 270 
K and 277 K. This strong temperature dependence and, 
particularly, the large DINS value reported for Kh at 271 
K (58% larger than at room temperature) are the main 
discrepancies found with our results for water. Our sim- 
ulations predict that under isothermal conditions Kh is 
larger in ice than in water while the values derived from 
DINS data at 269 K forecast that Kh is much lower in 
ice Ih than in water. 

In the face of the discrepancies found between our 
simulation results for Kh and those derived from neu- 
tron scattering experiments, we have related Kh to other 
thermodynamic properties pursuing to gain a deeper un- 
derstanding of the implication of such differences. We 
find that the maximum in the density of water is not 
related to any significant effect in the proton kinetic en- 
ergy at atmospheric pressure. In the analysis of C v our 
simulations display reasonable results for both water and 
ice Ih at 271 K. The Kh contribution to C p amounts to 
about 4 J mol _1 K _1 in water (5% of C p ). However, the 
temperature dependence of Kh as derived from DINS 
experiments leads to a dramatic negative contribution 
to C p of about -10 4 J mol -1 K -1 . This value can not 
be directly compared to experimental data, but both its 
sign and magnitude seem to be outside any reasonable 
physical boundary. At this point it is clear that further 
independent experiments and theoretical calculation are 
necessary to clarify this question. 

The calculated values of Kq and Kh in water and ice 
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at isothermal conditions show that the largest results are 
always associated to the solid phase. This fact determines 
the sign of the isotopic shifts in the melting temperature 
upon isotopic substitution of either H or O atoms. In 
both cases the model predictions are in agreement to ex- 
periment. The temperature dependence of Kh in ice Ih 
has shown satisfactory agreement to results derived from 
PIMD simulations based on ab initio DFT calculations 
using the SIESTA code ^21 

Summarizing, it is clear that an empirical potential 
model of water can not reproduce quantitatively all ther- 
modynamic properties of the real system. However, it is 
the combined analysis of several thermodynamic proper- 
ties such as the proton kinetic energy Kh, the intramolec- 
ular and intermolecular contributions to Kh, the heat 
capacity C p , the isotopic shifts in melting temperatures, 



and the density-temperature curve of water, what let us 
believe that the employed model provides realistic results 
for the temperature dependence of Kh in water and ice 
Ih at atmospheric pressure. 
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